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ABSTRACT 

The goal of this paper is to attempt to give some insight and guidelines on the 
application of nonlinear dynamic theory to the better understanding of steady-state 
numerical solutions and nonlinear instability in algorithm development for nonlinear 
differential equations that display genuinely nonlinear behavior in computational sci- 
ences and, in particular, computational fluid dynamics (CFD). This stems from the 
fact that, although the study of nonlinear dynamics and chaotic dynamics for nonlin- 
ear differential equations and for discrete maps have independently flourished rapidly 
for the last decade, there are very few investigators addressing the issue on the con- 
nection between the nonlinear dynamical behavior of the continuous systems and the 
corresponding discrete map resulting from finite-difference discretizations. This issue 
is especially vital for computational sciences since nonlinear differential equations in 
applied sciences can rarely be solved in closed form and it is often necessary to replace 
them by finite dimensional nonlinear discrete maps. In addition, it is also important to 
realize that these nonlinear discrete maps can exhibit a much richer range of dynamical 
behavior than their continuum counterparts. 

Furthermore, it is also very important to identify some of the implications of what 
happens when linear stability breaks down for problems with genuinely nonlinear behav- 
ior. Studies indicate that for relatively simple nonlinear ordinary differential equations 
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(ODEs) and well-known time- discretization with modest step-sizes some schemes can 
converge to a spurious (false) steady-state solution in a deceptively “smooth” manner. 
In some instances, spurious steady states may appear below the linearized stability 
limit of the schemes, and consequently computation may lead to erroneous results. Our 
preliminary studies on partial differential equations (PDEs) also show that much of 
nonlinear dynamic (e.g. chaotic) phenomena have a direct relation for problems con- 
taining nonlinear source terms such as the reaction-diffusion, the reaction-convection or 
the reaction-convection-diffusion equations. Here our object is neither to provide the- 
ory nor to illustrate with realistic examples the connection of the dynamical behavior 
of practical PDEs with their discretized counterparts, but rather to give insight into 
the nonlinear features unconventional to this type of study and to concentrate on the 
fundamental ideas. Thus, in order to bring out the special properties, the illustrations 
center on simple scalar differential equation (DE) examples in which the exact solutions 
of the DEs are known. 
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I. INTRODUCTION 


While the applied computational fluid dynamicists are busy developing numerical 
solver computer codes, grid generation codes, three-dimensional graphical stereo dis- 
plays, and stretching the limits of the faster supercomputers in the world to numerically 
simulate the various 3-D complex aerodynamic configurations [1], there is a group of 
applied mathematicians, physicists, chemists, biologists, applied mechanicians, and me- 
teorologists who are involved in a new science called “chaotic dynamics” (or nonlinear 
dynamics). The science of chaotic dynamics has cut across many traditional scientific 
disciplines for the last decade since chaotic dynamics is a science of the everyday world. 
It offers a way of seeing order and pattern where formerly only the random, the erratic, 
and the unpredictable were present. It explains much of the genuinely nonlinear phe- 
nomena that were once unexplainable. See references [2-10] for an introduction to this 
subject. 

Nonlinear Dynami cs & Chaotic Dynamics: Before the birth of chaotic dynamical the- 
ory, traditional study of nonlinear dynamics belonged to the applied mechanics disci- 
plines of mechanical engineering. Modern nonlinear dynamics (since the late seventies) 
includes chaotic dynamics. Thus, unless otherwise stated, the term nonlinear dynamics 
and chaotic dynamics are used interchangeably. That is, nonlinear dynamics includes 
chaotic dynamics and vice versa. 

Loosely speaking, the study of asymptotic behavior (steady-state solutions) of non- 
linear differential equations (DEs) and nonlinear discrete maps (difference equations) 
and how the asymptotes change as parameters of the system are varied is most of- 
ten referred to as nonlinear dynamic analysis and chaotic dynamic theory. Topics in 
this area include bifurcation theory, period doubling cascades resulting in chaos, etc. 
Stable chaotic solutions (chaotic attractors) may be defined loosely and simply as sta- 
ble asymptotes that have infinite period and yet are still bounded. It is emphasized 
here that unless otherwise stated, all DEs and discrete maps are nonlinear and consist 
of system parameters, and the terms discrete maps and difference equations axe used 
interchangeably. 

Types of Dynamical Systems: Consider an ordinary differential equation (ODE) and a 
partial differential equation (PDE) of the forms 

du 

— =aS{u), (1.1) 


du df(u) 
dt ^ dx 


d 2 v 


+ aS(u), 


(1.2) 


where a and e are parameters and 5 is a nonlinear function in u and is independent 
of q (and e). The function f(u) can be linear or nonlinear in u. An ODE of this 
form in which t does not appear explicitly in S is called an autonomous dynamical 
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system. One can also consider a function S which is nonlinear in v and depends on 
i. ODEs of this type are called nonautonomous dynamical systems and they are more 
difficult to analyze; see references [5,8] for a discussion. The analysis would be more 
complicated if S = 5(u, a) is nonlinear in both u and a. In this case, the DE is not only 
nonlinear in the dependent variable u (and independent variable t), but also nonlinear 
in the parameter space a. One can also consider systems that depend on more than one 
parameter and/or systems of equations of the above type. 

Next consider nonlinear discrete maps (nonlinear difference equations) of the forms 

U n+1 =u n + D(u n ,u n ~\r), (1.3) 

and 


«; + ’ = «? + G(„7,«7 ±1 ,r). (i.4) 

Here r is a parameter, and D is nonlinear in u n and v n ~ 1 and linear or nonlinear 
in the parameter space r. The situation is similar for the function G. One can also 
consider discrete systems that depend on more than one parameter. A typical example 
is a discrete map arising from a finite- difference approximation of DEs such as (1.1) or 
(1.2). For the ODE, the resulting discrete maps might be nonlinear in a as well as the 
time step At, depending on the ODE solvers. For the PDE, again depending on the 
differencing scheme, the resulting discretized counterparts can be nonlinear in a, At, 
the grid spacing Ax and the numerical dissipation parameters even though the DEs 
consist of only one parameter or none. 

One can also consider discrete maps (scalar or system) of the forms 


u 


n+J 


= u n + D(u 


n-f k 


.,U 


.71 — 1 


, Tj , 7*2, 


where are positive integers and rj,r 2 ,...,r m are parameters, and 


(1.5) 


u 


n41 


n , p/.n+fc n n-/ n-ffc 


u n 


,tt” ',ri, r 2 ,...,r m ). (1.6) 


Again, (1.6) can depend on more than the three indices j,j ± 1. Systems (1.4) and (1.6) 
are sometimes referred to as a partial-difference equation. The dynamical behavior of 
(1.4) and (1.6) can be many orders of magnitude more difficult than (1.3) and (1.5). 
Any of the systems ( 1 .1 )-( 1 .6 ) are examples of dynamical systems. 


Important Consideration : It is emphasized here that discrete maps, regardless of their 
origin, are dynamical systems on their own right. It is also important to distinguish the 
following five types of discrete maps: 

1. Discrete maps arise naturally in physical sciences. They commonly arise through 
the inability to measure populations at all points in space and time [5,10,11] in popula- 
tion dynamics. They can also arise through the study of periodic excitation of dyn ami cal 
systems [12,13] in applied mechanics. 
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2. Discrete maps arise from Poincare sections in ODEs [2]. 

3. Discrete maps arise from discrete approximations of ODEs. 

4. Discrete maps (partial-difference equations) arise from temporal and spatial finite 
difference approximations of PDEs. 

5. Discrete models arise from the “Inverse Problems of Nonlinear Dynamics” in time 
series analysis of observable data or experiments [9]. 

Discrete maps of types 1 and 5 sometimes might not have any relationship with a spe- 
cific continuum DE. As a matter of fact, there might be no concrete associated governing 
equations (continuum or otherwise) to start with for type 5 except the surrogated dis- 
crete map arising from the time series analysis. Type 2 arises naturally from the study 
of dynamical behavior of nonlinear ODEs. However, types 3 and 4 have an intimate 
link (but with a different tie than type 2) between the original governing continuum 
DEs and their discretized counterparts. Furthermore, it is important to distinguish the 
complexity involved in the analysis of types 3 and 4. Type 4 involves spatial as well as 
temporal dynamical behavior. 

Note that for disc rete maps of types 3 and 4, even though the DEs might be linear in 
the parameter space, depending on the numerical methods, the discretized counterparts 
might be linear or nonlinear in that parameter space. In addition, extra parameters 
which may be linear or nonlinear can also be introduced by the scheme as noted in 
the paragraph after equation (1.4). An important concept is that even though the 
DE does not depend on any parameter, its discretized counterpart does depend on at 
least one parameter. As can be seen in the subsequent sections, the nature of the 
dynamical behavior of these discrete maps is strongly influenced by properties of the 
numerical method and the types and forms of nonlinearity on the DEs. Furthermore, 
when dealing with nonlinear conservation laws of PDEs, the dynamical behavior of the 
discretized counterparts is also strongly influenced by elements such as conservation and 
nonlinearity of the schemes, and treatment of the source terms [14-18]. These issues are 
very crucial for the existence of spurious steady-state numerical solutions which will be 
explained in a later section. Here the term nonlinear scheme refers to a case where the 

resulting discrete maps are nonlinear when applied to scalar constant coefficient linear 
DEs. 

Objectives: Our ultimate objective is to conduct long term basic research on the inter- 
disciplinary field of integrating the theory of nonlinear dynamics with computational 
sciences and, in particular, with computational fluid dynamics (CFD). This new ap- 
proach to CFD is extremely difficult and complex to analyze. A summary of the diffi- 
culty involved was discussed in Yee [14] and will be elaborated in sections IV and V. Our 
immediate goal is to study the behavior of spurious steady-state numerical solutions for 
nonlinear DEs and the dynamical behavior of this type of numerical solutions. Even 
within this frame work, the subject is still very young, board, difficult and unfamiliar 
to computational scientists as well as researchers working in nonlinear dynamics and 
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nonlinear physics. 

The intent of this paper is to give a flavor of the subject, to familiarize the reader 
in computationl sciences with this new and exciting area, and most of all, to explain 
through simple illustrations why it is so important for computational scientists to learn 
about the subject. Some challenging topics for future research are also proposed. 

Because of the complexity involved, there is a vast difference in the degree of difficulty 
on the study of the subject between the discretized counterparts of nonlinear ODEs and 
the discretized counterparts of nonlinear PDEs. In order to achieve our final goal of 
studying the dynamical behavior of numerical methods for nonlinear PDEs that arise 
from, e.g., computational fluid dynamics (CFD), we have to first fully understand this 
subject on the time discretization and later link this knowledge to the study of both 
the temporal and spatial nonlinear dynamical behavior of finite-difference methods for 
nonlinear PDEs of the nonhomogeneous hyperbolic and parabolic types. 

Therefore, the content of this paper will concentrate on the dynamical behavior of 
time discretization for ODEs or systems of ODEs obtained from time- splitting [19] 
or method of lines [20] for PDEs, and emphasis will be placed on its implication for 
algorithm development in CFD and computational sciences in general. Hopefully this 
will be part I of a series of many future research papers to come under the same topic. 
Our companion paper [21] studies the dynamical behavior of the class of explicit Runge- 
Kutta methods in detail. The intent of this paper is to not only serve as a study 
of the state-of-the-art of nonlinear dynamical behavior of ODE solvers, but also more 
importantly to serve as an introduction and to present new results to motivate this vast, 
new yet unconventional concept. Thus the mission of this paper is not to provide the 
answer or theory or to illustrate the connection of dynamical behavior of practical PDEs 
to their discretized counterparts, but rather to gain insight into the nonlinear features 
unconventional to this type of study and concentrate on the fundamentals. In order 
to bring out the new features, the illustrations concentrate on the simple scalar DEs 
examples in which the exact solutions of the DEs are known. 

Outline: The outline of the paper is as follows: First, a brief background, motivation and 
basic ideas will be given. Then some typical characteristics of dynamical systems with 
genuinely nonlinear behavior will be discussed. Next, the dynamical behavior of discrete 
maps arising from time discretization of ODEs will be studied and the main results and 
their implications for computational sciences will be described. Studies on discrete maps 
arising from finite-difference approximations of PDEs will not be elaborated. Rather, 
the level of complexity involved and state-of-the-art study on this subject will be briefly 
described. The paper will conclude with a few recommendations. Remarks will be given 
on the popular misconception of residual test for convergence in steady-state solution via 
the “time-dependent” approach and the popular misconception of the use of the “Inverse 
Problems of Nonlinear Dynamics” to analyze the dynamical behavior of time series data 
from a computer code in an attempt to learn about the true physical solution behavior 
of the governing PDEs. This application of time series analysis can be misleading and 
a wrong conclusion can be reached if the practitioner does not know by other means 
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other than the numerical solutions the exact solution behavior of the PDEs . 
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II. MOTIVATION & RELEVANCE 


As discussed in the introduction, dynamical systems occur in the form of DEs and 
discrete maps. In order to motivate why the study of numerical analysis will not be 
complete without the utilization of the nonlinear dynamic approach, and to convey to 
practitioners in computational sciences the importance of distinguishing the difference 
between weakly nonlinear problems and genuinely nonlinear problems, this section is 
devoted to a discussion of the typical behavior of dynamical systems with genuinely non- 
linear behavior and the basic characteristic difference in dynamical behavior between 
DEs and discrete maps in general. This discussion leads to the key elements of this 
paper, namely: (1) to establish the connection between the DEs and their discretized 
counterparts and (2) to convey to computational scientists how one should change the 
traditional way of thinking and practices when dealing with genuinely nonlinear prob- 
lems. 


2.1. Typical Characteristics of Dynamical Systems 
with Genuinely Nonlinear Behavior 

The terms “nonlinear behavior” and “genuinely nonlinear behavior” are used quite 
often in the literature and there seems to be no unified exact definition or meaning 
[9]. Here these terms are used for nonlinear dynamical systems that exhibit mainly the 
following characteristics. 

(1) The study of nonlinear dynamics most often emphasizes the importance of ob- 
taining a global qualitative understanding of the character of the system’s dynamics 
since local analysis is not sufficient to give the global behavior of genuinely nonlinear 
dynamical systems. As a matter fact, this is one of the major reasons why sometimes 
it required orders of magnitude more work than solving their linear counterparts. 

(2) Unlike linear or weakly nonlinear problems, the solutions of genuinely nonlinear 
DEs and discrete maps are strongly dependent on initial data, boundary conditions and 
system parameters. 

(3) Only genuinely nonlinear dynamical systems can have chaotic behavior and one 
of the striking characteristics of chaotic behavior is sensitivity of the solution to initial 
data. This characteristic is independent of whether the dynamical system is a continuum 
or a discrete map. 

From here on, the terms “dynamical systems with genuinely nonlinear behavior” and 
“genuinely nonlinear dynamical systems” are used interchangeably. For convenience, 
the word “genuinely” is omitted in most parts of the paper. 
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2.2. Typical Difference in Dynamical Behavior of 
ODEs and Discrete Maps 

The study of discrete maps is the discrete analog to the study of ODEs, as the 
study of recursion formulas is a discrete analog to the study of series expansions of 
functions. Much of the theory of ODEs can carry over to discrete maps with some 
slight modifications. However, there are new phenomena occurring in discrete maps 
which are absent in differential systems [22,23,12,13]. 

With repect to the topographical behavior, there are new kinds of behavior of tra- 
jectories in the neighborhood of equilibrium points of discrete maps. The behavior of 
separatrices associated with a saddle type of equilibrium point for a nonlinear differ- 
ence system is far more complicated than the behavior of separatrices for a differential 
system. See Yee, Hsu and Hsu et al. [12,13,24,25] for details and examples. 

With respect to similar equation types, the minimum number of first-order nonlinear 
autonomous ODEs is three for the existence of chaotic phenomena. However, a simple 
scalar first-order difference equation [26-30] like the logistic map 

v n 

v n+1 = ftv n (l — ), (i a parameter, (2-1) 

or its piecewise linear approximation [31] 


v 


n+1 


= MV n , 
= 


v n < 1 
1 < t’ n < 3 


= fi(4 — v n ) 3 < v n . 


( 2 . 2 ) 


possesses very rich dynamical behavior such as period-doubling cascades resulting in 
chaos. Equation (2.2) has the same behavior as (2.1) except that simple closed form 
asymptotic solutions of all periods can be obtained. These characteristic trade differ- 
ences between ODEs and discrete maps are very general. The discrete maps can arise 
from any of the five types as discussed in the introduction. It is in this spirit that we 
say that discrete maps can exhibit a much richer range of dynamical behavior than 
DEs. The next two sections focus on the typical difference and connection between the 
dynamical behavior of ODEs and their discretized counterparts. 


2.3. Background and Motivation 

Spurious asymptotic numerical solutions such as chaos were observed by Ushike [32] 
and Brezzi et al. [33] on the leapfrog method for the logistic ODE 


du 

~dt 


— au( 1 — u ). 


(2.3) 
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In reference [34], Schreiber and Keller discussed the existence of spurious asymptotic 
numerical solutions for a driven cavity problem. Some related studies are reported in 
(35J. 


Spurious solutions of Burgers’ equation and channel flows have been studied and com- 
puted in [36-38]. Many other investigators in the computational sciences (e.g. [39-43]) 
have observed some kind of strange or chaotic behavior introduced by the numerical 
methods but were not able to explain systematically the source, the cause of their 
results, or most of all the implication and impact in practical applications in compu- 
tational sciences. Due to the popularity of searching for chaotic phenomena, it is very 
trendy to relate inaccuracy in numerical methods with the onset of chaos. It is em- 
phasized here that inaccuracy in long time integration of discrete maps resulting from 
finite discretization of nonlinear DEs comes in other forms prior to the onset of chaotic 
phenomena. Stable and unstable spurious steady states and spurious periodic numeri- 
cal solutions set in before chaotic behavior occurs. These spurious asymptotes of finite 
period are just as inaccurate as chaotic phenomena as far as numerical integration is 
concerned. In other words, the prelude to chaotic behavior is the key element that we 
want to stress (i.e., before the the onset of chaos or a divergent solution) since the result 
of operating the time step beyond the linearized stability limit is not always a divergent 
solution in genuinely nonlinear behavior; spurious steady-state solutions can occur. As 
can be seen at a later section, this behavior is more difficult to detect than chaotic 
phenomena in practical computations. 

Recently, it has been realized by numerical analysts that numerical methods for ODEs 
and PDEs can be considered as dynamical systems. Several papers [44,45] on numeri- 
cal methods as dynamical systems have appeared in recent years. These investigators 
studied the dynamical behavior of the different ODE solvers per se without relating its 
close tie with the ODEs themselves. Although the study of chaotic dynamics for non- 
linear differential equations and for discrete maps have independently flourished rapidly 
for the last decade, there are very few investigators addressing the issue of the con- 
nection between the nonlinear dynamical behavior of the continuous systems and the 
corresponding discrete map resulting from finite difference discretizations. This issue 
is especially vital for computational sciences since nonlinear differential equations in 
applied sciences can rarely be solved in closed form and it is often necessary to replace 
them by finite dimensional nonlinear discrete maps. Most often, typical applied scien- 
tists rely on numerical methods to give insight into the solution behavior of nonlinear 
DEs. It is not always clear how well a numerical solution can mimic the true physics of 
problems that possess genuinely nonlinear types of behavior. 

Why is there such a need to study the connection between the continuum and its 
discretized counterparts for CFD applications? This stems from the fact that current 
supercomputer power can perform numerical simulations on virtually any simple 3-D 
aerodynamic configuration and, due to the limited available experimental data, the 
applied engineers are relying on or trusting the numerical simulations whole heartedly 
to help design our next generation aircraft and spacecraft. However, many of these 
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applied scientists are still using linear analysis as their guide to study highly nonlinear 
equations, and most often they are not aware of the limitations and pitfalls of many of 
the numerical procedures. Furthermore, most of the numerical algorithms in use operate 
under the accuracy and stability limit guidelines of the linearized model equation. It 
is only appropriate to analyze nonlinear problems with the nonlinear approach; i.e., by 
the nonlinear dynamic approach. 

The unique dynamical property of the separate dependence of solutions on initial 
data for the individual nonlinear DE and its discretized counterpart is especially impor- 
tant for employing a “time-dependent” approach to the steady state with given initial 
data in hypersonic CFD. In many CFD computations, the steady-state equations are 
PDEs of the mixed type and a time-dependent approach to the steady state can avoid 
the complication of dealing with elliptic-parabolic or elliptic-hyperbolic types of PDEs. 
However, this time-dependent approach has created a new dimension of uncertainty. 
This uncertainty stems from the fact that in practical computations, the initial data 
are not known and a freestream condition or an intelligent guess for the initial condi- 
tions is used. In particular the controversy of the “existence of multiple steady-state 
solutions” through numerical experiments will not be exactly resolved until there is a 
better understanding of the separate dependence on initial data for both the PDEs and 
the discretized equations. 


2.4. Connection Between the Dynamical Behavior of the 
Continuum and Its Discretized Counterpart 


Aside from truncation error and machine round-off error, a more fundamental dis- 
tinction between the continuum and its discretized counterparts is new behavior in the 
form of spurious stable and unstable asymptotes created by the numerical methods. 
This is due to the fact that nonlinear discrete maps can exhibit a much richer range of 
dynamical behavior than their continuum counterparts as discussed in section 2.2. Some 
instructive examples will be given in section III. These new phenomena were partially 
explored by the University of Dundee group [46-54], Sanz Serna [55], Iserles [56,57] and 
Stuart [58-62]. Their main emphasis was on phenomena beyond the linearized stability 
limit. The main contribution of our current study is (1) the occurrence of spurious 
steady-state numerical solutions below the linearized stability limit of the scheme for 
genuinely nonlinear problems, (2) the strong dependence of numerical solutions on the 
initial data, as well as other system parameters of the DEs such as boundary conditions 
and numerical dissipations terms, and (3) the implications for practical computations 
in hypersonic CFD. 

Before discussing the numerical examples, the next two subsections will give an overall 
summary of our current findings (integrating with other relevant recent results). The 
discussion is divided into steady-state solutions and asymptotes of any period, and 
transient solutions. 
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2.4.1. Steady-state Solutions and Asympotes of Period Higher Than One: 


Table 2.1 shows the possible stable asymptotic solution behavior between DEs (ODEs 
or PDEs) and their discretized counterparts. Some of the phenomena will be supported 
by simple examples in section III. The main connection between the DEs and their 
discretized counterparts is that steady-state solutions of the continuum are solutions 
of the discretized counterparts but not the reverse. Their main difference is that new 
phenomena are introduced by the numerical methods in the form of spurious stable 
and unstable asymptotic solutions of any period. In the past, the phenomena of spu- 
rious asymptotes were observed largely beyond the linearized stability of the schemes. 
Some numerical analysts and applied computational scientists were not alarmed and 
were skeptical about these phenomena since, theoretically, one is always guided by the 
linearized stability limit of the scheme. However, this resaoning is only valid if one is 
solving a scalar nonlinear ODE and the initial data are known. Another important 
concept is that the result of operating with time steps beyond the linearized stability 
limit is not always a divergent solution; spurious steady-state solutions and spurious 
asymptotes of higher period can occur. 

Our current study also indicated that depending on the form of the nonlinear DEs, all 
ODE solvers can introduce spurious asymptotic solutions of some period or all periods. 
However, the most striking result is that for certain schemes and depending on the form 
of the nonlinear DEs, spurious steady states can occur below the linearized stability 
limit. See section III and our companion paper [21] for more details. 

Another important factor is that associated with the same (common) steady-state 
solution, the basin of attraction (domain of attraction) of the continuum might be 
vastly different from the discretized counterparts. This is due entirely to the separate 
dependence and sensitivity on initial and boundary conditions for the individual system. 
The situation is compounded by the existence of spurious steady states and asymptotes 
of period higher than one and possibly chaotic attractors. 

Here the basin of attraction of a dynamical system is the domain for which the set of 
initial conditions time asymptotically approaches a specific asymptote. Figures 2.1 and 
2.2 show the basins of attraction of two popular ODE dynamical systems. Figure 2.1 
shows the multiple stable steady states and their basins of attraction for the damped 
pendulum equation 


du 



(2.4a) 


— = — tv - sin(u) (2.4b) 

dt 

for t = 0.4. Figure 2.2 shows the multiple steady states and their basins of attraction 
for the simple predator-prey equation 

Y = — 3u 4- 4 u 2 — uv/2 — u 3 , (2.5a) 
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d v 

— — — 2.1V + uv. 
at 


(2.5b) 


where u is the population of the prey and v is the population of the predator. These 
figures are taken from Parker and Chua [8] and were generated by the use of a variable 
time step Runge-Kutta-Fehlberg method with built in accuracy check (if the numerical 
solutions are approximating the true solution of the ODE). See reference [8] for details. 
These figures, although generated numerically, with the built in accuracy check the fixed 
points and basins of attraction coincide with the ODEs. The stable fixed points of the 

damped pendulum equation are 2n7r, n = 0,1, The unstable fixed points (saddles) 

are (2n + l) 7 r. The separatrices of the saddle points divide the phase plane into the 
different basins of attraction for the corresponding stable fixed points. The fixed points 
of the predator and prey equation are slightly less regular than the damped pendulum 
equation. Figure 2.2 shows two saddle points at u = 1, v = 0 and ti = 3, v = 0, one 
stable focal point at u = 2.1, v = 2 and one stable nodal point at u = 0, v = 0. Again 
the separatrices of the saddle points divide the phase plane into the basins of attraction 
for the corresponding stable fixed points. 

Intuitively, in the presence of spurious asymptotes, the basin of the true steady states 
(steady states of the DEs) can be separated by the basins of attraction of the spurious 
asymptotes and interwoven by unstable asymptotes, whether due to the physics (i.e., 
present in both the DEs and the discretized counterparts) or spurious in nature (i.e., 
introduced by the numerical methods). 

For PDEs, another added dimension is that even with the same time discretization 
but different spatial discretizations or vice versa, the basins of attraction can also be 
extremely different. However, mapping out the basins of attraction for any nonlinear 
continuum dynamical system other than the very simple scalar equations relies on nu- 
merical methods. The type of nonlinear behavior and the dependence and sensitivity 
to initial conditions for both the PDEs and their discretized counterparts make the 
understanding of the true physics extremely difficult when numerical methods are the 
sole source. Under this situation, how can one delineate the numerical solutions that 
approximate the true physics from the numerical solutions that are spurious in nature? 
Hopefully, with our simple illustrations in section III, we can demonstrate the impor- 
tance of the current subject and, most of all, stress the importance of knowing the 
general dynamical behavior of asymptotes of the schemes for genuinely nonlinear scalar 
DEs before applying these schemes in practical calculations. 


2.4.2. Transient or Time-Accurate Solutions : 

It is a common misconception that inaccuracy in long time behavior poses no conse- 
quences on transient or time-accuate solutions. This is not the case when one is dealing 
with genuinely nonlinear DEs. For genuinely nonlinear problems, due to the possible 
existence of spurious solutions, larger numerical errors can be introduced by the nu- 
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merical methods than one can expect from local linearized analysis or weakly nonlinear 
behavior. The situation will get more intensified if the initial data of the DE is in the 
basin of attraction of a chaotic transient [63-65] of the discretized counterpart. This is 
due to the fact that existence of spurious asymptotes transact wrong behavior in finite 
time. In fact, it is possible the whole solution trajectory is likely to be erroneous. 

We’d like to end this section with a direct quote from Sanz- Serna and Vadillo’s paper 
[55]. This quote indicates the danger of relying on linearized stability and convergence 
theory in analyzing nonlinear dynamical problems. Reference [55] is one of the few pa- 
pers trying to convey to numerical analysts the flavor of the powerful “nonlinear dynamic 
approach”. Hopefully, with the current discussion, we can convey to computational fluid 
dynamicists the flavor of the importance of the “nonlinear dynamic approach” in CFD 
analysis. 

“Assume that the convergence of a numerical method has been established; it is 
still possible that for a given choice of At, or even for any such a choice, the qual- 
itative behaviour of the numerical sequence t/°, u 1 , u n , ... be competely differ- 
ent from that of the theoretical sequence u(t 0 ), u(t 3 ), ••• u (*n)i ••• This discrep- 
ancy which refers to n tending to oo, At fixed cannot be ruled out by the conver- 
gence requirement, as this involves a different limit process (namely At tending to 

0 ). 

. The fact that analyses based on linearization cannot ac- 
curately predict the qualitative behaviour of u n for fixed At should not be surprising: 
there is a host of nonlinear phenomena (chaos, bifurcations, limit cycles ...) which 
cannot possibly be mimicked by a linear model.” 
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III. THE ODE CONNECTION 


In this section, we review some of the fundamentals and available theory and discuss 
our maior results. The discussion will have some overlap with our companion paper 
[ 21 ]. 


3.1. Preliminaries 


Consider an autonomous nonlinear ODE of the form 

^ = qS(u), (3.1) 

where a is a parameter and S(u) is nonlinear in u. For simplicity of discussion, we 
consider only autonomous ODEs where a is linear in (3.1); i.e., a does not appear 
explicitly in S. 

A fixed point u * of an autonomous system (3.1) is a constant solution of (3.1); that 
is 


S(u*) = 0. (3.2) 

Note that the terms “equilibrium points”, “critical points”, “stationary points”, 
“asymptotic solutions” (exclude periodic solutions for the current definition), “steady- 
state solutions” and “fixed points” are sometimes used with slightly different meanings 
in the literature, e.g., in bifurcation theory. For the current discussion and for the ma- 
jority of nonlinear dynamic literature, these terms are used interchangeably. We might 
want to mention that certain researchers reserve the term “fixed point” for discrete 
maps only. 

Consider a nonlinear discrete map from finite discretization of (3.1) 

u n+I = u n + I)(u n ,r), (3.3) 

where r = aAt and D(v n .r) is linear or nonlinear in r depending on the ODE solvers. 
Here the analysis is similar if D is a nonlinear function of u n+p , p = 0, 1 , ..., m. Examples 
to illustrate the dependence on the numerical schemes for cases where D is linear or 
nonlinear in the parameter space will be given in the subsequent section. 

A fixed point u* of (3.3) (or fixed point of period 1) is defined by u n+1 = u n , or 

u* =u’ + D(u*,r) (3.4a) 


or 


D(u*,r) ---- 0. 


(3.4b) 
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One can also define a fixed point of period p, where p is a positive integer by requiring 
that u n + p = u n or 

u*=E p (u*,r ) but u*^E k (u*,r) for 0 < k < p. (3.5) 

Here, E p (u*,r ) means that we apply the difference operator E p times, where 
E(u n ,r) = u n + D(u n ,r). For example, a fixed point of period 2 means u n+2 = u n or 

u* = E(E(u\r)). (3.6) 

In this context, when dealing with discrete systems, the term “fixed point” without 
indicating the period means “fixed point of period 1” or the steady-state solution of 
(3.3). 

In order to illustrate the basic idea, the simplest form of the Ricatti ODE, i.e., the 
logistic ODE (2.3) with 


S(u) = u(l — it) 

is considered. For this ODE, the exact solution is 


(3.7) 


u(t) = 


u° 4- (1 — u°)e _OIt ’ 


(3.8) 


where u° is the initial condition. The fixed points of the logistic equation are roots of 
u*(l — tt*) = 0; it has two fixed points u* = 1 and u* = 0. 

To study the stability of these fixed points, we perturb the fixed point with a distur- 
bance and obtain the perturbed equation 


^ = q 5( u * +0- 

Next, S(u* -|- £) can be expanded in a Taylor series around «*, so that 


(3.9) 


dt 


— Q 


S(0 + s„(*0( + is..(Of , + 


(3.10) 


where S u (u') = | .• Stability can be detected by examining a small neighborhood 

of the fixed point provided if for given a, u* is not a hyperbolic point [3,7,9] (i.e., if the 
real part of aS u (u m ) 0). Under this condition £ can be assumed small, its successive 
powers £ 2 ,£ 3 , ... can normally be neglected and the following linear perturbed equation 
is obtained 


-±=aS u ( U ')Z- (3-11) 

The fixed point u m is asymptotically stable if aS u {u") < 0 whereas u* is unstable if 
aS u (u") > 0. If q5„(u*) = 0, a higher order perturbation is necessary. 
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If we perturb the logistic equation around the fixed point with a > 0, one can find 
that u* = 1 is stable and u* — 0 is unstable. It is well known that the general asymptotic 
solution behavior of the logistic ODE is that for any u° > 0, the solution will eventually 
tend to it* = 1. Figure 3.1 shows the solution behavior of the logistic ODE. 

Now, let us look at three of the well known ODE solvers. These are explicit Euler 
(Euler, forward Euler), leapfrog and Adam-Bashforth. For the ODE (3.1) with S(u ) = 
u(l — a), the dynamical behavior of their corresponding discrete maps is well established. 
The explicit Euler is given by 


u n+1 = u n + rS(u n ), (3.12) 

and it is after a linear transformation, the well known logistic map [26-30]. The leapfrog 
scheme can be written as 


u n+l = u n-l + 2 rS(u n ), 


(3.13) 


and it is a form of the Henon map [32]. The Adam-Bashforth method given by 


u 


n+J 



3S(u n ) - 5(u n_1 ) , 


(3.14) 


is again a variant of the Henon map and has been discussed by Pruffer [44] in detail. 

We can determine fixed points of the discrete maps (3.12)-(3.14) and their stability 
properties in a similar manner as for the ODE. It turns out that all three of the discrete 
maps have the same fixed points as the ODE (3.1) — a desired property which is im- 
portant for obtaining asymptotes of nonlinear DE numerically. Here we use asymptotes 
to mean fixed points of any period. 

The corresponding linear perturbed equation for the discrete map (3.3), found by 
substituting u n = v m + f n in (3.3) and ignoring terms higher than is 


c +1 = r [1 + AtD u {u*, At)] . (3.15) 

Here the parameter a of the ODE has been absorbed in the parameter At due to the 
assumption that a does not appear explicitly in S(u). For stability we require 


|1 + AtD v (u\At)\ < 1. (3.16) 

Again, for |1 + AtD u (u*, Af)| = 1, higher order perturbation is necessary. For a fixed 
point of period p the corresponding linear perturbed equation and stability criterion are 

= £"££(“', At). (3-17) 

and 

\E*{u*,At)\ < 1, (3.18a) 
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with 


££(u n ,A<) = ^-E{u n+ ^\At)...^-E{u n ,At) (3.18b) 

du du 

For S(u ) = ti(l — u), the stability of the stable fixed points of period 1 and 2 for 
discrete maps (3.12)-(3.14) with r = aAi are 

Explicit Euler: 

u* = 1 stable if 0 < r < 2 

period 2 stable if 2 < r < y/6. 

Leapfrog: 

u* = 1 unstable for all r > 0 

chaotic solution exist for all r no matter how small 

Adam-Bashforth: 

u* — 1 stable if 0 < r < 1 

period 2 stable if 1 < r < s/l. 

Figure 3.2 shows the stable fixed point diagram of period 1,2, 4, 8 by solving numer- 
ically the roots of (3.12) for S(u) = u(l — u). The r axis is divided into 1,000 equal 
intervals. The numeric labelling of the branches denotes their period. The subscript E 
on the period 1 branch indicates the stable fixed point of the DE. 

Two of these three examples serve to illustrate that the result of operating with a 
time step beyond the linearized stability limit of the stable fixed points of the nonlinear 
ODEs is not always a divergent solution; spurious asymptotes of higher period can 
occur. This is in contrast to the ODE solution, where only a single stable asymptotic 
value u* = 1 exists for any a > 0 and any initial data u° > 0. It is emphasized here that 
these spurious asymptotes, regardless of the period, stable or unstable, are solutions in 
their own right of the discrete maps resulting from a finite discretization of the ODE. 


3.2. Spurious Steady-State Numerical Solutions 


For the previous three ODE solvers, we purposely picked the type of schemes that do 
not exhibit spurious fixed points [56] but allow spurious fixed point of period higher than 
1. In this section, we discuss the existence of spurious steady-state numerical solutions. 
Again, it is emphasized here that these spurious steady states, stable or unstable, are 
solutions in their own right of the resulting discrete maps. Consider two second-order 
Runge-Kutta schemes, namely, the modified Euler (R-K 2) and the improved Euler (R- 
K 2), the fourth-order Runge-Kutta method (R-K 4), and the second and third-order 
predictor-corrector method [66-68] of the forms 
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Modified Euler (R-K 2) method: 

u n+1 = w n + rS(u n + ^rS n ^ S n = S(v n ) 

Improved Euler (R-K 2) method: 

u n+1 = u n + T - 1 5 [u n + S(u n + rS n )] | 

R-K 4 method: 

+ 2k2 + 2&3 + k 4 ^ 

ki = S n 

k 2 =s(u n +^rk 1 ) 
k 3 = S (u n + Kkij 

k 4 = S^u n + irA: 3 ^) (3.21) 

Predictor-corrector method of order m: 
u (0) = u n + rS n 

+ J 5 n + 5 (fc) |, * = 0,l,...,m-l 

u n+i = u n + - 5 n + 5 (Tn " 1) . (3.22) 

2 

Using the same procedures, one can obtain the fixed points for each of the above 
schemes (3.18) - (3.22). Figures 3.3 - 3.7 show the stable fixed point diagrams of period 
1,2,4 and 8 for these five schemes for S(u) = u(l — u). Some of the fixed points of lower 
period were obtained by closed form analytic solution and/or by a symbolic manipulator 
such as MAPLE [69] to check against the computed fixed point. The majority are 
computed numerically [2,8]. The stability of these fixed points was examined by checking 
the discretized form of the appropriate stability conditions. Again the axis is divided 
into 1,000 equal intervals. The numeric labelling of the branches denotes their period, 
although some labels for period 4 and 8 are omitted due to the size of the labelling areas. 
The subscript E on the main period one branch indicates the stable fixed point of the 
DE while the subscript S indicates the spurious fixed points introduced by the numerical 
scheme. Spurious fixed points of period higher than one axe obvious and are not labeled 



(3.18) 

(3.19) 
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except for special cases. Note that these diagrams, which appear in most parts as solid 
lines are actually points, which are only apparent in areas with high gradients. 

To contrast the results, similar stable fixed point diagrams are also computed for 
S(u) = u(l - u)(b - it), 0 < b < 1. See figures 3.8 - 3.14. The stable fixed point for 
the ODE in this case is u* = b and the unstable ones are u* = 0 and u* = 1. For any 
0 < u° < 1 and any a > 0, the solution will asymptotically approach the only stable 
asymptote of the ODE u* = b. 

Note that contrary to the DE, the maximum number of stable and unstable fixed 
points (real and complex) for each scheme varied between 4 to 16 for S(u) = u(l — u) 
and 9 to 81 for S(u) = u(l — u)(b — u), depending on the numerical methods and the r 
value. The domains of all of the fixed point diagrams are chosen so that they cover the 
most interesting part of the scheme and ODE combinations. Notice that asymptotes 
might occur in other parts of the domain as well. 

Aside from the striking difference in topography in the stable fixed point diagrams 
of the various methods and ODE combinations, all of these diagrams have one similar 
feature; i.e., they all exhibit spurious stable fixed points, as well as spurious stable 
fixed points of period higher than one. Although in the majority of cases, these occur 
for values of r above the linearized stability limit, this not always the case, as in the 
modified Euler scheme applied to the logistic ODE and du/dt = au(l — u)(b — u), 
0 < b < .5, and the R-K 4 applied to the logistic DE. For these two methods and ODE 
combinations, stable spurious fixed points occur below the linearzed stability limit. In 
some of the instances, these spurious fixed points are outside the interval of the stable 
and unstable fixed points of the ODEs. Others not only lie below the linearized stability 
limit but also in the region between the fixed points of the DEs and so could be very 
easily achieved in practice. 

One might argue that for the ODEs that we are considering, it is trivial to check 
whether an asymptote is spurious or not. For example, if a is a spurious asymptote of 
period one, then S(u) ^ 0. The main purpose of the current illustration is to set the 
baseline dynamical behavior of the scheme so that one can use it wisely in other more 
complicated settings such as when nonlinear PDEs are encountered in which the exact 
solutions are not known. Under this situation, spurious asymptotes could be computed 
and mistaken for the correct steady-state solutions. 

Note that fpr the modified Euler method, spurious fixed points of higher periods and 
chaotic attractors as well as spurious steady states occur below the linearized stability 
limit. Let fi be the basin of attraction of the fixed point of the ODE and let r* be 
the corresponding linearized stability limit value of the scheme. Then there exists a 
portion of the basin D denoted by Q c in which D c C D and an interval of r with 
0 < t < r* which actually belongs to the basin of attraction of the chaotic attractor of 
the discretized counterparts. There also exist some other D p C D and an interval of r 
with 0 < r < r* and p > 1 an integer, which actually belongs to the basin of attraction 
of a stable asymptote of period p of the corresponding discrete map. This leads to the 
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issue of the dependence of solutions on initial data which will be a subject of the next 
subsection. 


3.3. Strong Dependence of Solutions on Initial Data 


For simple nonlinear ODEs that we are considering, the fixed point diagram is ex- 
tremely useful for the understanding of the dynamics of the DEs and their discretized 
counterparts. However, when fixed points of higher periods and/or complex nonlinear 
equations are sought, searching for the roots and testing for stability of highly compli- 
cated nonlinear algebraic equations can be expensive and might lead to inaccuracy. 

Equally useful for understanding the dynamics are the bifurcation diagram and basin 
of attraction of fixed points for both the DEs and the difference schemes. The bifurcation 
diagram for the one-dimensional discrete maps displays the iterated solution u n vs. r 
after iterating the discrete map for a given number of iterations with a chosen initial 
condition (or multiple initial conditions) for each of the r parameter values. 

Bifurcation is broadly used to describe significant qualitative changes that occur in the 
orbit structure of a dynamical system as the system parameters are varied. In general, 
bifurcation theory can be divided into two general classes, namely, local and global. 
Local bifurcation theory is concerned with the bifurcation of fixed points of nonlinear 
equations and discrete maps. Global bifurcation studies phenomena away from the fixed 
points. It studies the interaction between different types of fixed points. One might 
define a bifurcation point as being any dynamical system which is structually unstable 
[3,8,9]. A fixed point is structurally stable if nearby solutions have qualitatively the 
same dynamics. The linearized stability limit of a fixed point of a scheme is the same as 
the bifurcation point in the corresponding bifurcation diagram of the resulting discrete 
map. 

For the numerical computations of the bifurcation diagrams with a given interval of r 
and a chosen initial condition (or multiple initial conditions), the r axis is divided into 
500 equal spaces. In each of the computations, the discrete maps were iterated with 
600 preiterations and the next 200 iterations were plotted for each of the 500 r values. 
The domains of the r and u n axes are chosen to coincide with the stable fixed point 
diagrams shown previously. For our current interest, it is not necessary to distinguish 
the difference between a stable fixed point of period 200 and a chaotic attractor. 

Figure 3.15 shows the bifurcation diagram of the Euler scheme applied to the logistic 
DE with an initial condition u° > 0. It is of interest to know that in this case the 
bifurcation diagram looks practically the same for any u° > 0. This is due to the fact 
that no spurious fixed points or spurious asymptotes of low periods exist for r < 2.627. 
Comparing the bifurcation diagram with figure 3.2, one can see that if we computed all 
of the fixed points of period up to 200 for figure 3.2, the resulting fixed point diagram 
would look the same as the corresponding bifurcation diagram (assuming 800 iterations 


21 



of the logistic map are sufficient to obtain the converged stable asymptotes of period 
upto 200 and a proper set of initial data are chosen to cover the basins of all of the 
periods in question). The numeric labelling of the branches in the bifurcation diagram 
denote their period, with only the essential ones labelled for identification purposes. 

In order to interpret the bifurcation diagram for other ODE and scheme combina- 
tions, some knowledge of the fixed point diagram is necessary, at least for the lower 
order periods. Otherwise, one cannot identify the exact periodicity of the asymptotes 
easily. As can be seen later, a “full” bifurcation diagram cannot be obtained efficiently 
without the aid of the stable and unstable fixed point diagram for schemes that exhibit 
spurious fixed points of any period, especially lower periods. In most cases, the un- 
stable asymptotes divide the domain into the proper basins of attraction for the stable 
asymptotes (spurious or otherwise), and at least one initial data point is used from each 
of the basins of attraction before a full bifurcation diagram can be obtained. 

In all of the fixed point diagrams 3.3 - 3.14, the bifurcation phenomena can be divided 
into three kinds. For the first kind, the paths (spurious or otherwise) resemble period 
doubling bifurcations (flip bifurcation) [2-5] similar to the logistic map. See figures 3.2, 
3.6 and 3.8 for examples. The second kind occurs, most often, at the main branch 1 £, 
with the spurious paths branching from the correct fixed point as it reaches the linearized 
stability limit, and quite often even bifurcating more than once (pitchfork bifurcation 
or supercritical bifurcation [70,7]), as r increases still further before the onset of period 
doubling bifurcations. See figures 3.4, 3.7, 3.9 - 3.11 and 3.13 for examples. The 
third kind again occurs most often at the main branch 1#. The spurious paths near the 
linearized stability limit of 1# would experience a transcritical bifurcation [3,7,9,70]. See 
figures 3.3, 3.5, 3.7 and 3.14 for examples. Notice that the occurrence of transcritical 
and supercritical bifurcations are not limited to the main branch l#. See figures 3.11 
- 3.14 for examples. The other commonly occurring bifurcation phenomenon is the 
subcritical bifurcation which was not observed in our two chosen S(u) functions. With 
a slight change in the form of our cubic function 5(u), a subcritical bifurcation can be 
achieved [70,3,7,9]. The consequence of the latter three bifurcation behaviors is that 
bifurcation diagrams calculated from a single initial condition u° will appear to have 
missing sections of spurious branches, or even seem to jump between branches. This is 
entirely due to the existence of spurious asymptotes of some period or more than one 
period, and its dependence on the initial data. This occurs even for the Euler scheme 
as depicted in figure 3.8. See section 3.4 for further discussion of these four types of 
bifurcation phenomena. 

Figures 3.16 - 3.18 show the bifurcation diagram by the modified Euler method for the 
logistic ODE with three different starting initial conditions. In contrast to the explict 
Euler method, none of these diagrams look alike. One can see the influence and the 
strong dependence of the asymptotic solutions on the initial data. Figure 3.19 shows 
the corresponding “full” bifurcation diagram, their earlier stages resembling the fixed 
point diagram 3.3. Figures 3.20 - 3.22 illustrate similar bifurcation behavior for the 
corresponding R-K 4 method. Figure 3.12 serves as an example to illustrate that the 
effect of overplotting a number of initial data, but not the appropriate ones, would not 
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be sufficient to cover all of the essential spurious branches. Figures 3.23 - 3.25 show a 
similar illustration for S(u) = u(l — u)(b — u), 0 < b < .5 by the improved Euler, R-K 
4 and the modified Euler method. The strong dependence of solutions on initial data 
is evident from the various examples in which this type of behavior is very common for 
genuinely nonlinear problems. 

In order to compute a “full” bifurcation diagram, we must overplot a number of 
diagrams obtained by the guide of the stable and unstable fixed point diagram as an 
appropriate set of starting initial data. In the case where the fixed point diagrams are 
extremely difficult to compute, a brute force method of simply dividing the domain of 
interest of the u n axis into equal increments and using these u n values as initial data 
is employed. The “full” bifurcation diagram is obtained by simply overplotting all of 
these individual diagrams on one. 

For completeness, figures 3.26 - 3.38 show the “full” bifurcation diagrams for the 
corresponding fixed point diagrams shown previously. Figures 3.36 and 3.37 show a 
blow up section of figures 3.34 and 3.35. Notice that the exact values of the initial 
data are immaterial as long as these values cover all of the basins of attraction of the 
essential lower order periods (i.e., at least one initial data point is used from each of 
the basins). Here, we use the term “full” bifurcation diagram to mean just that. No 
attempt has been made to compute the true full bifurcation diagram since this is very 
costly and involves a complete picture of the basins of attraction for the domain of 
interest in question. 


3.4. Classification of ODE solvers 
(According to the Existence of Spurious Fixed Points) 


In reference [56], lserles studied the stability of ODE solvers for nonlinear autonomous 
ODE via the dynamical approach. He proved that linear multistep methods (LMM) [66- 
68] that give bounded values at infinity always produce correct asymptotic behavior, but 
it is not the case with Runge-Kutta methods and some predictor-corrector methods. He 
demonstrated that the Runge-Kutta and predictor- corrector methods may lead to false 
asymptotes. However, he did not discuss the possibility of these spurious asymptotes 
existing below the linearized stability limit. 

For implicit LMM, he assumed the resulting nonlinear algebraic equations are solved 
exactly. He also showed the influence of nonlinear algebraic solvers on the size of stability 
regions for implicit LMM. His conclusion was that the standard nonlinear algebraic 
solver — the modified Newton- Raphson method 
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can drastically degrade the region of stability limit as compared to the Newton-Raphson 
method 
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On the other hand, the direct iteration method 
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converges only if the step sizeis of the same order of magnitude as that required for an 
explicit method. Thus the advantage of using an implicit method to enhance stability 
is lost. Here for clarity of notation, when iteration procedures are involved, u n is used 
in place of u n of the previous section. 


The implications of behavior detailed in Iserles’ work [56] range far beyond pure ODE. 
For most CFD application, the use of implicit time discretization to “time” march the 
solution to steady state is very common. The resulting nonlinear algebraic systems 
are solved by either noniterative linearization [71,14] or by some kind of iterative or 
relaxation procedures. Very often, applied computational fluid dynamicists experience 
a non-convergent solution where the residual will decrease only so far before reaching 
a plateau with a time step larger than the explicit method. Therefore the behavior 
observed in Iserles’ work could explain the degradation in the stability of the implicit 
scheme in practice. Indeed, even though the mechanisms involved are far more compli- 
cated than those studied here, elements such as spatial discretization dynamical behavior 
and nonlinear coupling effect for systems, could well be an explanation. 

More recently, Iserles and Sanz- Serna [57] established conditions for using a variable 
step size analysis to avoid spurious fixed points in a class of Runge-Kutta methods. 

Looking at the problem from another perspective, it is very useful to find the cause 
of the existence of spurious asymptotes by looking at the form and properties of the 
resulting discrete maps, regardless of the methods. We have the following two observa- 
tions. 


(1) Assume that the only parameter that was introduced by a numerical method 
is At. Then from Iserles results and our current investigation, one obvious necessary 
condition for the existence of spurious steady states of ODE solvers for (3.1) is the intro- 
duction of nonlinearity in the parameter space At. This is evident from our examples 
and general analysis. For example, if At (or r) is linear in (3.3), then (3.3) can be 
written as 


u n+1 = u n + cr5'(u n ), c a constant of the scheme. (3.26) 

Therefore any fixed point of (3.3) is a fixed point of (3.1). Without lost of generality, a 
similar proof applies to the resulting difference operator D from a p time level scheme. 

(2) The second observation is that one can classify the types of spurious steady 
state in the form of bifurcation theory near a bifurcation point or a bifurcation limit 
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point. Figures 3.39 and 3.40 show the definition of the various types of branching 
points and the stability of solution in the neighborhood of branch points. In other 
words, the classification is according to the onset of spurious asymptotes of subcritical, 
supercritical or transcritical bifurcations. See figure 3.41 for the definition of the three 
types of phenomena. 

Assume an ODE solver introduces nonlinearity in the parameter space At for (3.1). 
Then a necessary and sufficient condition for the occurrence of spurious steady states 
below the linearized stability limit on the main branch 1# (stable fixed points of the 
DE) is that a transcritical or subcritical bifurcation of the types shown in figures 3.42 
and 3.43 exist at the bifurcation point or near a bifurcation limit point. It is emphasized 
here that the existence of spurious fixed points of higher period can be independent of 
the existence of spurious steady states (fixed points of period 1). 

A detailed analytical analysis on the existence of transcritical, subcritical and super- 
critical bifurcations for the class of Runge-Kutta methods can be found in our com- 
panion paper [21]. Figures 3.44 - 3.54 illustrate the onset of different types of spurious 
steady states by showing the stable and unstable fixed points of periods 1 and 2, and 
the types of bifurcation phenomena for the modified Euler, Improved Euler and R- 
K 4 and the predictor-corrector schemes of order 2 and 3 for S(u) = u(l — u) and 
S(u ) = u(l — u)(b — u), 0 < b < .5. In order to illustrate the different behavior in 

an uncluster fashion, not all of the periods 1, 2 and branching points are labeled. It is 
interesting to see the manner in which the onset of the different types of bifurcations 
occur, in particular, the birth of the different types of bifurcations away from the 1 £ 
branches. 


3.5. Basins of Attraction 


Due to the separate dependence and sensitivity on initial data for the individual 
DEs and the discretized counterparts, in conjunction with the existence of spurious 
steady states and asymptotes of higher periods, even associated with the same (common) 
steady-state solution, the basin of attraction of the continuum might be vastly different 
from the discretized counterparts. 

Take for example, S(u) = u(l — u). The only stable fixed point of the logistic ODE 
is v = 1. The entire domain of the real u n -axis is divided into two basins of attraction 
for the ODE independent of any positive a. Now if one numerically integrates the ODE 
by the modified Euler method, extra stable and unstable fixed points can be introduced 
by the scheme depending on the value of r. That is for certain ranges of the r values, 
the u n -axis is divided into four basins of attraction. But of course for other ranges of r, 
higher period spurious numerical solutions exist, more basins of attraction are created 
within the same u^-axis range, etc. Stable and unstable fixed point diagrams such as 
figures 3.44 - 3.54 are very useful in the division of the u”-axis into different basins of 
lower periods. 
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3.6. Systems of ODEs 


As can be seen from the previous sections, the rich and complicated dynamical be- 
havior of discrete maps resulting from finite discretization of simple nonlinear scalar 
autonomous ODEs is very enlightening, educational and useful in giving some indica- 
tions of the strange behavior encountered in practice. One would naturally ask how 
highly coupled nonlinear first-order autonomous systems complicate the issue. After 
all, these types of systems occur naturally in physical science and engineering fields. 
Examples are 

(1) second or higher order nonlinear scalar autonomous or nonautonomous ODEs 
arising from mechanical systems, 

(2) meteorology, 

(3) chemical reaction equations arising from chemistry, 

(4) system of ODEs arising from the method of lines approach in reaction- diffusion, 
reaction-convection and reaction-convection-diffusion equations. 

Future work will be directed towards investigation into the nonlinear dynamical effect 
of using ODE solvers for nonlinear system of ODEs. Here, we do not attempt to give 
a detailed discussion on this subject, but rather indicate some of the implications from 
our experience as well as what is availiable in the literature. 

First, the coupling of first-order nonlinear systems arising from a higher-order scalar 
nonlinear ODE is very different from the truly nonlinear coupling on systems of first- 
order ODEs. This difference carries over to their discretized couterparts. Second, due 
to the nonlinear coupling effect, whatever is observed in the nonlinear scalar case will 
definitely exists in the coupled system case in a more complex manner. Even with the 
help of the center manifold theorem [2-5], nonlinear systems of higher than three first- 
order ODEs are still extremely difficult to analyze. One major factor in analyzing the 
associated discrete maps from finite discretization of the continuum is that when three 
or more time levels of ODE solvers are used, even though the continuum is a first-order 
scalar autonomous ODE, the resulting discrete maps are (p — 1 )th-order, where p is the 
time level. One can extrapolate the complexity involved if nonlinear coupled systems of 
higher-order ODEs were discretized by p-time levels of ODEs solvers. Some aspects and 
implications of numerical integration of second- and third-order ODEs are discussed in 
references [39,40,72]. Some of our preliminary numerical experiments agree with the 
above general conclusion. 


3.7. Suitability to the Type of Computational Environment 


The main approach that we use in this paper is to establish the necessary mathemati- 
cal reasoning and then to support this reasoning with extensive numerical experiments. 
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Our current study on one-dimensional scalar nonlinear dynamical equations which con- 
sist of a single parameter indicates that the understanding of the nonlinear effects en- 
countered when applying finite-difference schemes to nonlinear differential equations is 
greatly aided by the analysis of bifurcation diagrams which record the values of succes- 
sive iterations for a range of parameters. Equally useful are diagrams showing the basins 
of attraction of equilibria, both those of the differential equations and the spurious at- 
tractors generated by the difference scheme. The generation of such diagrams, however, 
is computationally expensive, especially for the basins of attraction where each point 
on the diagram represents a different choice of parameters for which many iterations of 
the scheme must be performed to determine its significance. 

In all the bifurcation diagrams, the computations were performed on the VMS VAX 
in double precision. Take for example, figure 3.34. Each dot on the plot represents a 
solution obtained by integrating the discretized equation 800 times with each of the 20 
prescribed initial data and each of the 500 equally space values of r = a At. In other 
words, we are integrating the same equation for 10,00Q different values of r and initial 
data combinations, and also iterating the same equation for each of these combinations 
with 800 iterations. The task can therefore be greatly enhanced by parallel computation, 
since essentially the same process needs to be applied to each point in a fine two- or 
three-dimensional array, each element representing a pixel on a high resolution screen or 
plotter. It is therefore a task highly suited to machines, such as the Connection machine, 
which have large numbers of processors enabling the entire region or subregions of the 
problem to be analysed in one pass rather than in a sequential point-by-point approach. 
The intensity of (repetitive) computing involved is too great to gain major benefit from 
machines such as the CRAY. 

For multidimensional systems consisting of several parameters, we envision that the 
intensity of repetitive computing to obtain a bifurcation diagram or a basin of attraction 
cannot be realized if it is not performed on a massively parallel computer such as the 
Connection machine. 
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IV. LEVEL OF COMPLEXITY FOR PDEs 


In order to systematically approach the subject of studying spurious steady-state 
numerical solutions of nonlinear nonhomogeneous hyperbolic and parabolic PDEs via 
the nonlinear dynamic approach, we propose to pursue the subject in three stages. First, 
we will attempt to obtain a full understanding of the subject for time discretization of 
ODEs. The investigation can give insight into numerical methods employing the Strang 
type of operator splittings or methods of lines approach for nonhomogeneous hyperbolic 
and parabolic PDEs. The second stage will involve the study of the discrete travelling 
wave solutions of the reaction-convection and reaction-convection-diffusion equations. 
The third stage will involve the study of the complete temporal-spatial discretizations 
of the reaction-convection and reaction-convection-diffusion equations. The last stage of 
the proposed plan is extremely difficult to analyze. Some aspects of full discretizations 
and discrete travelling wave solutions were investigated by [46-54, 58-62,73,74,10]. 

The question now is in what specific area will this approach advance the state-of-the- 
art in CFD. Our preliminary study indicated that many existing results for nonlinear 
dyn ami cal systems such as chaos, bifurcations, and limit cycles (closed periodic or- 
bits [5]) have a direct application to problems containing nonlinear source terms such 
as the reaction- diffusion, reaction-convection or the reaction-convection-diffusion equa- 
tions. Also they have a direct application to most of the nonlinear shock-capturing 
methods such as the total variation diminishing (TVD) schemes [14,75-78]. With the 
advent of increasing demand for numerical accuracy, stability, efficiency, and uniqueness 
of numerical solutions in modeling such equations, an interdisciplinary approach for the 
analysis of these systems and schemes is needed. Besides it is a common practice in CFD 
to employ a time-dependent approach to achieve steady state. The separate dependence 
of solutions on initial data and system parameters for the individual PDE and its finite- 
difference equations is the crucial element in determining how well a numerical solution 
can mimic the true physics of the problem. 

The following is an attempt to give a flavor of the subject and at the same time provide 
a justification for the importance of this subject area in CFD algorithm development 
for our next generation aerodynamics needs. 

4.1. Model Equations 

One of the recent areas of emphasis in CFD has been the development of appropri- 
ate finite-difference methods for nonequilibrium gas dynamics in the hypersonic range 
[14,78-81]. A nonlinear scalar reaction- diffusion model equation would be of the form 


du _ d 2 u 
dt £ dx 2 

a nonlinear scalar reaction-convection model equation would be of the form 


-f aS(u), e,a system parameters, 


(4.1) 
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= ct5(u), 


du df(u) 
dt + dx 


(4.2) 


and a nonlinear scalar reaction-convection-diffusion model equation would be of the 
form 


du dfju) _ 
dt dx £ dx 2 


+ aS(u). 


(4.3) 


Here f(u) a linear or nonlinear function of u. The nonlinear source term (or the reaction 
term) 5(u) can be very stiff. Note that phenomena such as chaos, bifurcations and limit 
cycles only relate to source terms .S'(u) which are nonlinear in u. Equation (4.3) can be 
viewed as a model equation in combustion or as one of the species continuity equations in 
nonequilibrium flows (except in this case, the source term is coupled with other species 
mass fractions). 


The above model equations are good starting points in the investigation of correlation 
between the theory of chaotic dynamical systems and uniqueness, stability, accuracy and 
convergence rate of finite- difference methods for CFD. 


4.2. Level of Complexity 


The main interest is to investigate what types of new phenomena arise from the 
numerical methods that are not present in the original nonlinear PDE, as a function 
of the stiff coefficient a, the diffusion coefficient e, and the time step At with a fixed 
(or variable) grid spacing Ax. The time step can vary greatly depending on whether 
the time discretization is explicit or implicit. More precisely, one wants to weed out all 
undesirable phenomena due to the numerical method (e.g., additional equilibrium points 
introduced by the time as well as spatial discretizations, degradation of the domain of 
attraction, etc.) and to identify whether the numerical method really describes the true 
solution of the PDE under prescribed initial and boundary conditions with a, e, the 
time step At and the grid spacing Ax being parameters. The study can be divided into 
steady and unsteady behavior with or without shock waves. 

The major stumbling block is that combustion-related and high speed hypersonic flow 
problems usually contain multiple equilibrium states and shock waves that are inherent 
in the governing equations. Furthermore, spurious equilibrium states can be introduced 
by the time differencing and/or the spatial differencing. In many instances the stable 
and unstable equilibrium states, whether due to the physics or spurious in nature, are 
interwoven over the domain of interest and are usually very sensitive to the initial 
conditions and the time steps (even when the chosen time step is within the linearized 
stability limit as indicated in our study) as well as variation of parameters such as angle 
of attack, Reynolds number and coefficients of physical and numerical dissipations and 
physical and numerical boundary conditions. 
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The sensitivity of numerical solutions to coefficients of physical and numerical dissipa- 
tions is evident from the study of Mitchell and Bruch on the reaction-diffusion equation. 
Their main result is that diffusion, which is usually perceived as having a stabilizing 
effect, is able to produce chaotic as well as divergent numerical solutions. Another in- 
teresting result due to Mitchell and Bruch was the production of chaos by decreasing 
the space increment or increasing the time increment. They showed that the addition of 
diffusion poses severe problems unless waves of constant speed c are assumed, in which 
case it reverts to an ODE with x + cl as the independent variable. The sensitivity 
of numerical solutions to numerical boundary condition procedures was discussed in 
[82,83]. 

On the subject of sensitivity and dependence of solutions on initial data, the basin 
of attraction might be very different between the PDE and the discretized counterpart. 
The basin of attraction might contract or be very different from the basin of attraction 
for the original PDEs depending on the numerical methods. In many instances, even 
with the same spatial discretization but different time discretizations, the basins of 
attraction can also be extremely different. One can extrapolate the complexity involved 
when the influence of the various temporal as well as spatial discretizations are sought 
on the basins of attractivity. 

Table 4.1 summarizes the level of complexity for a systematic approach to these types 
of PDE. The check mark on each type of PDE and approach indicate the ones where 
some work has been done on this subject. The majority are credited to the University 
of Dundee group [46-54] and some related theory by A. Stuart [58-62]. 


4.3. Involvement in the Study of Full Discretization of PDE 


Consider a three-level explicit time differencing and a three-point spatial differencing 
of the reaction-convection-diffusion equation (4.3) of the form 


it 


n+ 1 


— it” + H » Wj , u j+i, u j-i i u j > u j+ 1 Ax), 


(4.4) 


where u” is the numerical solution at t = n£t and x = j Ax. Then the study of the 
asymptotes of (4.4) amounts to the study of fixed point behavior of period p in time and 
period q in space, denoted by (p, q), where p and q are integers. Here the fixed point of 
the partial-difference equation (4.4) is defined in a slightly more complicated way than 
for the ODE. 


For example, a fixed point of period (1,1) is defined as = u" and a fixed point 

of period (2,1) is defined as u"+ 2 = u”. However, a fixed point of period (1,2) is defined 
as ~ U J’ Thus, in general a fixed point of period (p,q) is defined as = u”. 

One can see that for p, q > 3, solving the resulting nonlinear algebraic equation is very 
involved, especially when physical boundary conditions and physical dissipation terms 
as well as numerical boundary conditions [82,83,34] and numerical dissipation [47] are 
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additional dimensions of consideration. Current available work involved the studies 
beyond the linearized stability limit of the schemes, and assumed the nonexistence of 
spurious fixed points of period (1,1). See references [46-54] for details. 


4.4. Influence in Dynamical Behavior by Property of the PDEs 
and Schemes, and Treatment of the Source Terms 


Although the general study of the dynamical behavior of parti ah difference equations 
for the conservation law [84,85] of (4.3) is an enormous task, if we can isolate certain 
restricted subsets of the PDEs and schemes in hand which are immune to the type of 
phenomena discussed in section III for time discretization as well as spatial discretiza- 
tion, then we can concentrate on the rest of the unknowns. 

As can be seen in section III, the nature of the dynamical behavior of the discretized 
counterparts is strongly influenced by properties of the numerical method and the types 
and form of nonlinear DEs. Here we want to study the influence on the dynamical be- 
havior of elements such as conservation and nonlinearity of the schemes, and treatment 
of the source terms [14-17,78-81] when nonlinear conservation laws of PDEs are sought. 

First, take the convection equation (4.2) with S(u ) = 0 and consider a conservative 
explicit scheme [76,14] which is consistent with the conservation law of the form 


u 


n + l 


n \ 

— u 3 — A 


hr.. t - h n i 

j+j j-\ 


(4.5) 


where A = At j Ax and h n , x are the numerical flux functions. For a two-time level and 

J ± 2 

five-point spatial scheme, x = h(tt^, 

We also can consider a two-parameter family of scheme 


1 1 \e 

< + ; 

; 1 + 0 ? 


l n +l 


h]±l 

J 2 


= u j ~ 


A(l-fl) 

1 4- t*> 


h\ x - h n 1 

J+2 }~7 


+ TTZ^- u T'y < 4 - 6 > 

where 0 < 6 < 1. When 0 = 0, the scheme is explicit and when 6 = u? + 1/2, the scheme 
is temporally second-order accurate. One can obtain (4.5) from (4.6) by setting 0 = 0 
and u) = 0. The time differencing belongs to the class of LMM. Under the assumption 
that this scheme is conservative and consistent with the conservation law, discrete map 
(4.6) will have no spurious steady-state numerical solution since consistency means 


=- f(u*). (4.7) 

Thus any steady-state solutions of (4.6) are steady-state solutions of the original PDE. 
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Now the situation is different when S'(u) 7^ 0. Under this situation, even if the 
same time and spatial discretization are employed, one still has to evaluate S properly. 
Here 5 is the function 5 evaluated at some proper average state u [14-17] for the 
full discretization that is consistent with the scheme [18], and achieves conservation at 
jumps. For a discussion on this subject, see references [78, 15-17] for details. The other 
crucial aspect is that when 5(u) ^ 0, a full investigation on the dynamical behavior 
of the temporal and spatial discretization is necessary. The knowledge gained from the 
finite-difference methods analysis for S(u) = 0 does not carry over to the S(u) ^ 0 case. 


4.5. Discrete travelling Waves 


Analysis of the dynamical behavior of the full discretization of nonlinear nonhomo- 
geneous PDEs of the hyperbolic and parabolic types is very involved. In this section, 
we look at a more restricted class of solutions — the discrete travelling wave solutions. 

Consider a reaction-diffusion equation 


du 

!k 


& + *<->• 


(4.8) 


Solution u(x,1 ) depend on the space variable x and on the time t. Every zero of 5(a) 
constitutes an equilibrium of the PDE. Then a travelling wave solution is a profile U (x) 
that travels along the x-axis with propagation speed A. Neither the shape of the wave 
nor the speed of propagation changes. To find travelling waves, we seek solutions 


resulting in an ODE 


u(x,<) = U(x — At), 


(4.9) 


U" + XU' + S(f/) = 0. (4.10) 

By solving this ODE, one can calculate asymptotic states for the PDE. Let U\ and 
U 2 be roots of S(u) and hence equilibrium solutions for both the PDE and ODE. The 
asymptotic behavior of solution U for x — ► ±00 determines the type of travelling wave. 
Every solution with 


lf(oo) = ttj (4.11a) 

U(-oo) = u 3 , (4.11b) 

with Uj / u 2 , is a front wave of the ODE. This corresponds to a heteroclinic orbit 
[3] of the ODE, connecting the two stationary points uj and u 2. Here for a second- 
order autonomous ODE (4.10), when distinct saddles are connected, one encounters 
a heteroclinic orbit; also a heteroclinic orbit may also join a saddle to a node or vice 
versa. Another type of special orbit is a homoclinic orbit. A homoclinic orbit connects a 
saddle point to itself and such orbits have an infinite period. Several heteroclinic orbits 
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may form a closed path called a homoclinic cycle. Both the heteroclinic and homoclinic 
orbits are of great interest in applications because they form the profiles of travelling 
wave solutions of many reaction-diffusion problems. See references [3,10,73,74] for a 
discussion. 

Similarly, one can study discrete travelling wave solutions for the finite discretization 
of (4.8). See references [73,74] for a discussion. Understanding of the discrete traveling 
wave solutions of the corresponding PDEs only gives insight into a very small subset of 
the dynamics of the PDEs. In most cases, it provides no information at all for the fully 
discretized equation. 
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V. IMPLICATIONS & RECOMMENDATIONS 


Due to the complexity of the large increase in system dimension and the involvement 
of multiple floating parameters for finite difference methods in PDEs, we are not certain 
that a similar systematic general result can be arrived at for more complex nonlinear 
systems. The main indication at this point is from our time discretization study. 


5.1. Results Drawn from the ODE Connection Study 

Our study illustrates a few very important implications which are very fundamental 
j n explaining what happens when linear stability breaks down for truly nonlinear prob- 
lems; i.e., equations that display genuinely nonlinear types of behavior. The important 
points are as follows: 

(1) There is sensitivity to initial data and strong dependence on discretization pa- 
rameters such as the time step and the grid spacing Ax. Dependence of solutions on 
initial condition is important for employing a time-dependent approach to the steady- 
state with a given initial condition and boundary conditions in hypersonic or combustion 
flows, especially when initial data of the governing PDE are not known. 

(2) Associated with the same (common) steady-state solution the basin of attraction 
of the DEs might be vastly different from the discretized counterparts. This is mainly 
due to the dependence and sensitivity on initial conditions and boundary conditions 
for the individual systems. In the absence of the influence of the initial and boundary 
conditions, the difference in the basins of attraction between the continuum and its 
discretized counterparts occurs even when an implicit LMM type of method is used 
unless the resulting nonlinear algebraic equations are solved exactly. 

(3) Nonunique steady-state solutions can be introduced by the spatial discretization 
even though the original PDEs might possess only an unique steady-state solution and a 
LMM type of time discretization is used so that no spurious steady-state exists in time. 
The tie between temporal and spatial dynamical behavior is more severe when one is 
dealing with the nonseparable temporal and spatial finite-difference discretization such 
as the Lax-Wendroff type, where the time and spatial difference cannot be separated 
from each other. The situation would be more complicated if the governing nonlinear 
PDE possesses more than one steady-state solution as well as the spurious ones that 
are purely due to the numerical method. 

(4) For certain time discretizations, spurious steady-state solutions may occur below 
the linearized stability limit- of the scheme. 

(5) The result of operating with a time step beyond the linearized stability limit is 
not always a divergent solution; spurious steady-state solutions can occur. 
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(6) There is a misconception that computational instability or inaccuracy can often 
be cured simply by making A t smaller. Other elements such as (1) - (5) above as well as 
the variation of the grid spacings, numerical dissipation terms and system parameters 
other than the time steps can interfere w T ith the dynamical behavior. 

(7) When linearized stability limits are used as a guide for a time step constraint 
for highly coupled nonlinear system problems, this time step might exceed the actual 
linearized stability limit of the coupled equations. Therefore all of the situations in 

(1) - (6) can occur. In particular, when one tries to stretch the maximum limit of the 
linearized allowable time step for highly coupled systems, most likely all of the different 
type of spurious branches of supercritical, subcritical and trancritical bifurcations can 
be achieved in practice depending on the initial conditions. That is why the occurrence 
of spurious steady-state solutions beyond the linearized stability limit is not just sec- 
ondary but might be as important as the occurrence of spurious steady states below the 
linearized stability limit. 


5.2. Recommendations 

It is of utmost importance to know the nonlinear dynamical behavior of the various 
schemes before their actual use for practical applications. Otherwise, it might be very 
difficult to asses the accuracy (spurious or otherwise) of the solution when the numerical 
method is the sole source of the understanding of the physical solutions. When in doubt, 
it is always safer to use schemes that do not produce spurious steady-state solutions for 
the nonlinear scalar case. Some examples of methods of this type in time discretization 
can be listed: 

(1) LMM [56] ODE solvers such as the explicit, implicit Euler, three-point backward 
differentation, etc. can be used. 

(2) One can use the “Regular” Runge-Kutta methods [57]. 

(3) Solving the nonlinear algebraic systems arising from implicit LMM method ex- 
actly would avoid spurious steady state numerical solutions. Otherwise, the type of 
iteration method in solving nonlinear algebraic systems can degrade the basin of attrac- 
tivity of implicit LMM [57]. 

The insight gained from time discretization will only give an indication in separable 
schemes or method of lines approaches. Also, the commonly used residual test [86-88] 
in the time-dependent approach to the steady state might be misleading. This is the 
direct consequence of what was indicated in section 5.1. The popular misconception of 
using the inverse problem of nonlinear dynamics to analyze a time series data from a 
finite difference method computer code in an attempt to learn about the true physical 
solution behavior of the continuum governing PDEs without knowing by other means 
the exact solution behavior of the PDEs other than the numerical solutions can also be 
misleading. These will be discussed in the next two sections. 
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5.3. Residual Test 


Consider a quasilinear PDE of the form 

du 

0 ^ == (5*1) 

where G is nonlinear in u , u z and u xx and a and e are system parameters. For simplicity, 
consider a two time level and a (p + g) point grid stencil of the form 


u 


»+i _ 


= ~ H(u] +q ,...,u^,...,u]_ p ,a,t,At,Ax) 


(5.2) 


for the PDE (5.1). Let U* , a vector representing (u* +? , ...,Uj_ p ) be a steady-state 

numerical solution of (5.2). It is a common practice in CFD to use a time dependent 
approach such as (5.2) to solve the steady-state equation G(u, u x ,u xx , a, e) = 0. The 
iteration is stopped when the residual H or some L 2 norm of the dependent variable u 
between two successive iterates is less than a pre-selected level. 


Aside from the various standard numerical error such as truncation error, machine 
round-off error, etc. [89], there is a more fundamental question on the validity of the 
residua] test and/or L 2 norm test. If the scheme happens to produce spurious steady- 
state numerical solutions, these spurious solutions would still satisfy the residual and L 2 
norm tests in a deceptively smooth manner. Moreover, aside from the spurious solutions 
issue, depending on the combination of time as well as spatial discretizations, it is not 
easy to check whether G(u* ,u*,u* x , a, e) — ► 0 even though H(U*,a,e,Ai,Ax) — * 0. 
This is contrary to the ODE case, where if u* is spurious in (1.1) then 5(u*) ^ 0. 
Among other factors, this is one of the contributing factors in the increase in magnitude 
of difficulty for analyzing the dynamical behavior of numerical methods for hyperbolic 
and parabolic PDEs. 


One might argue that one can judge the accuracy of the scheme by comparing the 
numerical solutions with more than one numerical methods and by doing a sequence of 
grid refinement and time step reductions. The latter approach might not be feasible at 
an acceptable cost. The former might not be foolproof if one does not know the dynam- 
ical behavior of the finite difference schemes being used. One important contributing 
factor on the use of the Laix-Wendroff types of schemes [90,91] is that these schemes 
are more accurate and sometimes more stable when operated on or near the linearized 
stability limit. 


5.4. The Inverse Problems of Nonlinear Dynamics 


The use of the inverse problem of nonlinear dynamics to analyze the dynamical behav- 
ior of time series data arising from experimental or observable data has received much 
attention in nonlinear physics as well as in many of the engineering disciplines. The 
approach is very useful for gaining some insights into the nonlinear dynamical behavior 
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in problems where experimental or observable data are the main source of informa- 
tion. Often the associated governing equations (continuum or otherwise) do not exist 
to start with. There has been an explosion of theory, numerical procedures and com- 
puter software addressing this rapidly growing direction [92-95]. There also has been 
much recent interest in forecasting algorithms that attempt to analyze a time series 
by fitting nonlinear models. The attractive feature of this approach is that when used 
correctly on the correct problems one can reduce the complexity of the problem from 
un-manageable higher dimensions to a very low dimension. It is therefore a natural 
tendency for practioners in computational sciences to apply this approach to analyze 
the dynamical behavior of time series data from a finite difference method computer 
code in an attempt to learn about the true physical solution behavior of the governing 
PDEs. This application of time series analysis can be misleading and can lead to a 
wrong conclusion if the practitioner does not know by other means the exact solution 
behavior of the PDEs other than from the numerical solutions. Examples of the use of 
this type of approach in CFD computations have been presented in references [96-98]. 
It can be seen from our study that the conclusions drawn from this type of time series 
analysis provide very little information, but rather can actually mislead one as to the 
true physics of the problem. 


VI. CONCLUDING REMARKS 


Spurious stable as well as unstable steady-state numerical solutions, spurious asymp- 
totic numerical solutions of higher period, and even stable chaotic behavior can occur 
when finite-difference methods are used to solve nonlinear DEs numerically. The oc- 
currence of spurious asymptotes is independent of whether the DE posseses a unique 
steady state or has additional periodic solutions and/or exhibits chaotic phenomena. 
The form of the nonlinear DEs and the type of numerical schemes are the determining 
factor. In addition, the occurrence of spurious steady states is not restricted to the time 
steps that are beyond the linearized stability limit of the scheme. In many instances, 
it can occur below the linearized stability limit. Therefore, it is essential for practi- 
tioners in computational sciences to be knowledgeable about the dynamical behavior of 
finite-difference methods for nonlinear scalar DEs before the actual application of these 
methods to practical computations. It is also important to change the traditional way 
of thinking and practices when dealing with genuinely nonlinear problems. 

In the past, spurious asymptotes were observed in numerical computations but tended 
to be ignored because they all were assumed to lie beyond the linearized stability lim- 
its of the time step parameter At. As can be seen from our study, bifurcations to and 
from spurious asymptotic solutions and transitions to computational instability not only 
are highly scheme dependent and problem dependent, but also initial data and bound- 
ary condition dependent, and not limited to time steps that are beyond the linearized 
stability limit. 

The symbiotic relation among all of these various factors makes this topic fascinating 
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and yet extremely complex. The main fundamental conclusion is that, in the absence of 
truncation and machine round-off errors, there are qualitative features of the nonlinear 
DE which cannot be adequately represented by the finite-difference methods and vice 
versa. The major feature is that convergence in practical calculations involved fixed At 
as n — + oo rather than At — + 0 as n — » oo. It should be emphasized that the resulting 
discrete maps from finite discretizations can exhibit a much richer range of dynami- 
cal behavior than their continuum counterparts. A typical feature is the existence of 
spurious numerical asymptotes that can interfere with stability, accuracy and basins of 
attraction of the true physics of the continuum. 
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Possible stable asymptotic solution behavior for DEs and their discretized 
counterparts. 

Phase portrait and basins of attraction of the damped pendulum equation 
(this figure is taken from reference [8]). 

Phase portrait and basins of attraction of the predator-prey equation (this 
figure is taken from reference [8]). 

Asymptotic solution behavior of the logistic ODE du/dt = ttt/.(l — u) for 
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Stable fixed points of periods 1,2, 4, 8 of the explicit Euler scheme for the 
logistic ODE du/dt = cm(l — u). 

Stable fixed points of periods 1,2, 4, 8 of the modified Euler (R-K 2) scheme 
for the logistic ODE du/dt = au(l — u). 

Stable fixed points of periods 1,2, 4, 8 of the improved Euler (R-K 2) scheme 
for the logistic ODE du/dt = au(l — u). 

Stable fixed points of periods 1,2, 4, 8 of the Runge-Kutta 4th-order (R-K 4) 
scheme for the logistic ODE du/dt = qu(1 — u). 

Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 2 for the logistic ODE du/dt = au(l — u). 

Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 3 for the logistic ODE du/dt = au(l — u). 

Stable fixed points of periods 1,2, 4, 8 of the explicit Euler scheme for the 
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Stable fixed points of periods 1,2, 4, 8 of the modified Euler (R-K 2) scheme 
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Stable fixed points of periods 1,2, 4, 8 of the improved Euler (R-K 2) scheme 
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Stable fixed points of periods 1,2, 4, 8 of the Runge-Kutta 4th-order (R-K 4) 
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Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 2 for the ODE du/di = au(l — u)(0.5 — u). 

Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 3 for the ODE du/dt = au(l - •u)(0.5 — u). 

Stable fixed points of periods 1,2, 4.8 of the modified Euler (R-K 2) scheme 
for the ODE du/dt = au(l - u)(b - u), b = 0.1, 0.2, 0.3, 0.4. 

Bifurcation diagram of the explicit Euler scheme for the logistic ODE du/dt = 
au( 1 — u). 

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE du/dt = au(l — u) with = 2.7. 

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE du/dt = au(l — u) with u° — 1.5. 

Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE du/dt = au(l — u) with u° = 0.25. 

“Full” bifurcation diagram of the modified Euler (R-K 2) scheme for the 
logistic ODE du/dt = au(l — u). 

Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for the 
logistic ODE du/dt = au( 1 — u) with u° = 0.5. 

Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for the 
logistic ODE du/dt = au(l — u) with multiple initial data. 

“Full” bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for 
the logislic ODE du/dt = att(l — u). 

Bifurcation diagrams of the improved Euler (R-K 2) scheme for the ODE 
du/dt = ou(l — u)(0.5 — u) for four different sets of initial input data. 

Bifurcation diagrams of the Runge-Kutta 4th-order (R-K 4) scheme for the 
ODE du/dt = au(l — u)(0.5 — u) for four different sets of initial input data. 

Bifurcation diagrams of the modified Euler (R-K 2) scheme for the ODE 
du/dt = «u(l — u)(0.4 — tt) for four different sets of initial input data. 

“Full” bifurcation diagram of the improved Euler (R-K 2) scheme for the 
logistic ODE du/dt = azt(l — u). 
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Fig. 3.27 “Fiill” bifurcation diagram of the Adam-Bashforth scheme for the logistic 
ODE du/di — au(l — u). 

Fig. 3.28 “Full” bifurcation diagram of the predictor-corrector scheme of order 2 for 
the logistic ODE dujdt = au(l — u). 

Fig. 3.29 “Full” bifurcation diagram of the predictor-corrector scheme of order 3 for 
the logistic ODE du/di = au(l — u). 

Fig. 3.30 “Full” bifurcation diagrams of the modified Euler (R-K 2) scheme for the 
ODE dujdt = au(l - u){b - u), b = 0.1, 0.2, 0.3, 0.4. 

Fig. 3.31 “Full” bifurcation diagrams of the explict Euler scheme for the ODE dujdt = 
au(l — u)(0.5 — u). 

Fig. 3.32 “Full” bifurcation diagram of the modified Euler (R-K 2) scheme for the 
ODE dujdt = au(l — tt)(0.5 — u). 

Fig. 3.33 “Full” bifurcation diagram of the Adam-Bashforth scheme for the ODE 
dujdt = cm(l — tt)(0.5 — u). 

Fig. 3.34 “Full” bifurcation diagram of the improved Euler (R-K 2) scheme for the 
ODE dujdt = qu(1 — u)(0.5 — u). 

Fig. 3.35 “Full” bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for 
the ODE dujdt = au(l — tt)(0.5 — u). 

Fig. 3.36 “Full” bifurcation diagram of the improved Euler (R-K 2) scheme for the 
ODE dujdt = au(l - u)(0.5 — ti) (enlarged). 

Fig. 3.37 “Full” bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for 
the ODE dujdt = ait(l — u)(0.5 — u) (enlarged). 

Fig. 3.38 “Full” bifurcation diagram of the predictor- corrector scheme of order 2 for 
the ODE dujdt = au(l — u)(0.5 — u). 

Fig. 3.39 Types of branching points. 

Fig. 3.40 Stability of solutions in the neighborhood of branch points, one-dimensional 

case. — stable, unstable a,b,c,d: limit (regular turning) point; 

e,f,g,h: bifurcation (double) points; ij,k,l: bifurcation-limit (singular turn- 
ing) points; m,n,o,p,q: additional possible cases when the dimension of u is 
greater than one (this figure is taken from reference [7]). 


48 


Fig. 3.41 Stability of steady-state solutions arising through three types of bifurcation 
phenomena ( — stable, - - - unstable). 

Fig. 3.42 Spurious fixed points arising from transcritical bifurcations. 

Fig. 3.43 Spurious fixed points arising from subcritical bifurcation. 

Fig. 3.44 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K 
2) scheme for the logistic ODE du/dt = au(l — u). 

Fig. 3.45 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K 
2) scheme for the logistic ODE du/dt — a«(l — u ). 

Fig. 3.46 Stable and unstable fixed points of periods 1,2 of the Runge-Kutta 4th-order 
(R-K 4) scheme for the logistic ODE du/dt = au(l - u). 

Fig. 3.47 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 
scheme of order 2 for the logistic ODE du/dt = au(l - u). 

Fig. 3.48 Stable and unstable fixed points of periods 1,2 of the predictor- corrector 
scheme of order 3 for the logistic ODE du/dt = au(l - u). 

Fig. 3.49 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K 
2) scheme for the ODE du/dt - mt(l - u)(0.5 - u). 

Fig. 3.50 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K 
2) scheme for the ODE du/dt = au(l - u)(0.5 - u). 

Fig. 3.51 Stable and unstable fixed points of periods 1,2 of the Runge-Kutta 4th-order 
(R-K 4) scheme for the ODE du/dt -- cra(l - u)(0.5 - u). 

Fi§" 3.52 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 

scheme of order 2 for the ODE du/dt = au (l - u)(0.5 - u). 

Fig. 3.53 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 

scheme of order 3 for the ODE du/dt = ou(l - tt)(0.5 - u). 

Fig. 3.54 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K 
2) scheme for the ODE du/dt = au(l - u)(0.2 - u). 

Table 4.1 Systematic approach - level of complexity. 
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SOLUTION 

TYPE 

ODEs or PDEs 

DISCRETIZED COUNTERPARTS 

# OF 

SINGLE 

SINGLE 

ASYMPTOTES 

SINGLE 

MULTIPLE 

OR 

STEADY-STATE 

MULTIPLE 

SAME # OF MULTIPLE 

SOLUTIONS 

MULTIPLE 

ADDITIONAL # OF MULTIPLE 

PERIODIC 

NO 

YES 

SOLUTIONS 

YES 

YES e- EXTRA) 

CHAOS 

NO 

YES 

YES 

YES ff EXTRA) 


Table 2.1 Possible stable asymptotic solution behavior for DEs and their discretized 
counterparts. 
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I. ODE CONNECTION: GAIN INSIGHT INTO TIME DISCRETIZATION OF PDEs 

SCALAR 

SYSTEM — time splitting or method of lines 

II. DISCRETE TRAVELLING WAVE: iy + r. £ ^ + 0 S(u) 

=== === ^ ==== at ax ax 2 

REACTION-DIFFUSION * 

SCALAR: REACTION -CONVECTION 

REACTION-CONVECTION-DIFFUSION 

III. FULL DISCRETIZATION (TEMPORAL AND SPATIAL): 


Table 4.1 Systematic approach - level of complexity. 


SCALAR: 

(S*» 


SCALAR: r 
(S#0) ( 


LINEAR SCHEME FOR SPATIAL DISCRETIZATION 
NONLINEAR SCHEME FOR SPATIAL DISCRETIZATION 
NONLINEAR SCHEME FOR SPATIAL DISCRETIZATION 
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Fig. 3.2 Stable fixed points of periods 1,2, 4, 8 of the explicit Euler scheme for the 
logistic ODE dujdt = au(l - u). 
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joints of periods 1,2, 4, 8 of the improved E 
~ic ODE du/dt = au(l - u). 
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Fig. 3.6 Stable fixed points of periods 1,2, 4, 8 of the predictor- corrector scheme of 
order 2 for the logistic ODE du/dt ~ au(l - u). 
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Fig. 3.9 Stable fixed points of periods 1,2, 4, 8 of the modified Euler (R-K 2) scheme 
for the ODE du/dt — att(l - u)(0.5 - u). 
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r = aAt 

Fig. 3.10 Stable fixed points of periods 1,2, 4, 8 of the improved Euler (R-K 2) scheme 
for the ODE du/di = au(l — u)(0.5 - u). 
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Fig. 3.11 Stable fixed points of periods 1,2, 4, 8 of the Runge-Kutta 4th-order (R-K 4) 
scheme for the ODE du/dt = au(l - u)(0.5 - u). 
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Fig. 3.12 Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 2 for the ODE du/dt = au(l - tt)(0.5 — u). 


65 



Fig. 3.13 Stable fixed points of periods 1,2, 4, 8 of the predictor-corrector scheme of 
order 3 for the ODE du/dt = au(l — u)(0.5 — u). 
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Fig. 3.16 Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE du/dt — au(l — u) with u° = 2.7. 
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Fig. 3.17 Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE du/di = au(l — u) with u° = 1.5. 
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Fig. 3.18 


Bifurcation diagram of the modified Euler (R-K 2) scheme for the logistic 
ODE dv/dt = ou(l — it) with u° = 0.25. 
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Fig. 3.20 Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for the 
logistic ODE du/dt = au(l — u) with u° = 0.5. 
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Fig. 3.21 Bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for the 
logistic ODE du/dt = au(l — u ) with multiple initial data. 
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Fig. 3.22 “Full” bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for 
the logistic ODE du/dt = au(l — u). 
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Fig. 3.23 


Bifurcation diagrams of the improved Euler (R-K 2) scheme for the ODE 
du/di = au(l - ti)(0.5 - u) for four different sets of initial input data. 
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Fig. 3.25 Bifurcation diagrams of the modified Euler (R-K 2) scheme for the ODE 
du/di — qu( 1 — u)(0.4 — u) for four different sets of initial input data. 
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Fig. 3.27 “Full” bifurcation diagram of the Adam-Bashforth scheme for the logistic 
ODE du/dkt = au(l — u). 
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Fig. 3.29 “Full” bifurcation diagram of the predictor-corrector scheme of order 3 for 
the logistic ODE du/dt = au(l — u ). 
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Fig. 3.30 “Full” bifurcation diagrams of the modified Euler (R-K 2) scheme for the 
ODE du/dt = qu( 1 - u)(6 - u), b = 0.1, 0.2, 0.3, 0.4. 





Fig. 3.31 “Full” bifurcation diagrams of the explict Euler scheme for the ODE du jdt = 
au(l — u)(0.5 — u). 
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ODE du/dt = aix(l — u)(0.5 — u). 
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Fig. 3.33 “Full’' bifurcation diagram of the Adam-Bashforth scheme for the ODE 
du/di = ctu( 1 — u)(0.5 — tt). 
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Fig. 3.35 “Full” bifurcation diagram of the Runge-Kutta 4th-order (R-K 4) scheme for 
the ODE du/dt = au(l — u)(0.5 — u). 


88 


1 . 0-1 



o H — i 1 1 1 1 1 1 

13.5 14.5 15.5 16.5 

r = a At 


Fig. 3.36 “Full” bifurcation diagram of the improved Euler (R-K 2) scheme for the 
ODE du/di = qu(1 — u)(0.5 — u) (enlarged). 
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Fig. 3.39 Types of branching points. 
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Fig. 3.40 Stability of solutions in the neighborhood of branch points, one-dimensional 
case. — stable, - - - unstable a,b,c,d: limit (regular turning) point; 

e,f,g,h: bifurcation (double) points; i,j,k,l: bifurcation-limit (singular turn- 
ing) points; m,n,o,p,q: additional possible cases when the dimension of u is 
greater than one (this figure is taken from reference [7]). 




Fig- 3.41 Stability of steady-state solutions arising through three types of bifurcation 
phenomena ( — stable, - - - unstable). 
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Stable and unstable 



. 3.45 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K 
2) scheme for the logistic ODE du/di — au(l — u). 
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. 3.46 Stable and unstable fixed points of periods 1,2 of the Runge-Kutta 4th-order 
(R-K 4) scheme for the logistic ODE du/di — au( 1 — u ). 






Fig. 3.48 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 
scheme of order 3 for the logistic ODE du/di = au(l - u). 
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Stable and unstable 
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Fig. 3.50 Stable and unstable fixed points of periods 1,2 of the improved Euler (R-K 
2) scheme for the ODE du/dt = au(l - u)(0.5 — u). 
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. 3.52 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 
scheme of order 2 for the ODE du/di = au(l — u)(0.5 - it). 
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Fig. 3.53 Stable and unstable fixed points of periods 1,2 of the predictor-corrector 
scheme of order 3 for the ODE du/di = au(l — u)(0.5 - u). 
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Fig. 3.54 Stable and unstable fixed points of periods 1,2 of the modified Euler (R-K 
2) scheme for the ODE du/dt = au(l - tt)(0.2 - u). 
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